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ABSTRACT 

Pulsar timing arrays (PTAs) can be used to search for very low frequency (10“®- 
10“^ Hz) gravitational waves (GWs). In this paper we present a general method for 
the detection and localization of single-source GWs using PTAs. We demonstrate the 
effectiveness of this new method for three types of signals: monochromatic waves as 
expected from individual supermassive binary black holes in circular orbits, GWs from 
eccentric binaries and GW bursts. We also test its implementation in realistic data 
sets that include effects such as uneven sampling and heterogeneous data spans and 
measurement precision. It is shown that our method, which works in the frequency 
domain, performs as well as published time-domain methods. In particular, we find 
it equivalent to the J>“Statistic for monochromatic waves. We also discuss the con¬ 
struction of null streams - data streams that have null response to GWs, and the 
prospect of using null streams as a consistency check in the case of detected GW 
signals. Finally, we present sensitivities to individual supermassive binary black holes 
in eccentric orbits. We find that a monochromatic search that is designed for circular 
binaries can efficiently detect eccentric binaries with both high and low eccentricities, 
while a harmonic summing technique provides greater sensitivities only for binaries 
with moderate eccentricities. 
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1 INTRODUCTION 


The direct detection of gravitational waves (GWs) will have profound implications in physics and astronomy. This may 
become possible within this decade. In the audio band (10-1000 Hz), second-generation km-scale laser interferometers, such as 
Advanced LIGO | |Harry|2010p , Advanced Virgo ( [Acernese et al.|2015p and KAGRA | |Somiya|2012[ ), are about to start scientific 
observations as early as the second half of 2015 and are expected to detect dozens of compact binary coalescence events per 
year when they achieve their design sensitivities in around 2020 (Abadie et al. 2010[ Aasi et al. 20131. In the nanohertz 
frequency range (1-100 nHz), high-precision timing observations of millisecond pulsars (pulsar timing arrays - PTAs) provide 
a unique means of detecting GWs. With the concept proposed decades ago ( Sazhin|[l978 Detweil^|1979 Hellings & Downs 


1983 Foster fc Backer]|1990 1, there are now three major PTA projects around the globe, namely, the Parkes Pulsar Timing 
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Array (PPTA; Manchester et aL]|2013 HobbspOlS I, the European Pulsar Timing Array ( Kramer &: Champion|[2013 ), and 
NANOGrav (McLaughlin 20131. While these PTAs have individually collected high-quality data spanning > 5 yrs for ~20 

Shannon et alj|20^l, they have also been combined to 


pulsars and produced some astrophysically interesting results (e.g., 


Manchester 


20131 aiming at significantly enhanced 


form the International Pulsar Timing Array (IPTA; Hobbs et al.||20i0 
sensitivities. 

The primary sources in the PTA frequency band are inspiralling supermassive binary black holes (SMBBHs). It is widely 
considered that a stochastic GW background due to the combined emission from a large number of individual SMBBHs over 
cosmological volume (e.g., Sesana|2013b Ravi et al.|2014 for recent work) provides the most promising target; indeed, some 
studies suggested that a detection of this type could occur as early as 2016 (Siemens et al. 20131. Analyses of actual PTA 


data have previously focused on a search for such a GW background, leading to more and more stringent constraints on the 
fractional energy density of the background (Jenet et al.|2006 Yardley et al.|20TT van Haasteren et al.|2011 Demorest et al. 


2013 Shannon et al. 20131. Over the past few years interest has grown substantially regarding the prospects of detecting 


single-source GWs using PTAs, for example, for individual SMBBHs (|Sesana et al. 


2012 


et al. 


2009 


2014|), for GW bursts ( Pitkin|[2012 I and for unanticipated sources (Gutler et al.|2014|. In the meantime, many data 


Lee et al. 


2011 


Ravi et ^|2015 1, for GW memory effects associated with SMBBH mergers (Setoj20M Cordes'fc'jenet|2oT2 


Mingarelli et al. 


Madison 


analysis methods have been proposed in the context of PTAs for single-source GW detection, for example, for monochromatic 
GWs emitted by SMBBHs in circular orbits ( Yardley et al.pOlO Babak&^^Sesana|20^ Ellis et al. 2012 EUisj|20^ Taylor 


|et al.|2014| |^^ng^^et^^alj20M Zhmrtalj201^, for GW memory effects ( vairHaastereiTfc Levin]201o| ' Wang^ al.|2oi5 1, and 

GW bursts ( EhunTToinmenpoTo Deng[2oT4| ). A few papers have presented searches in real PTA data for GWs from circular 
SMBBHs and yielded steadily improved upper limits on the GW strain amplitude ( Yardley et al.pOlO Arzoumanian et al. 
2014 Zhu et aL]|2014l. While circular binaries emit GWs at the second harmonic of the orbital frequency, eccentric binaries 


radiate at multiple harmonics. Jenet et al. (20041 first derived the expression of pulsar timing signals produced by eccentric 


binaries and developed a framework in which pulsar timing observations can be used to constrain properties of SMBBHs. 

In our previous work ( Zhu et ^[2014 1, we performed an all-sky search for GWs from SMBBHs in circular orbits using 
the latest PPTA data set. Here we adapt the network analysis method used in the context of ground-based interferometers 
(e.g., Pai et al. 2001 Wen & Schutz 2005 Wen & Schutz 2012 Wen 2008 Klimenko et al. 2008 Sutton et al. 20101 as a 


general method for detection and localization of single-source GWs using PTAs. In particular, we consider the following types 
of sources: (1) SMBBHs in circular orbits; (2) eccentric SMBBHs; and (3) GW bursts. We demonstrate the effectiveness of 
this method using synthetic data sets that contain both idealized and realistic observations. 

The organization of this paper is as follows. In sectionj^we briefly review the features of pulsar timing signals caused by 
single-source GWs and discuss the signal models used in this work. In section we present the mathematical framework of 
our method and propose practical detection statistics. The relation to the method used in our previous work is also discussed. 
Using idealized simulations we show examples of detection, localization and waveform estimation in section]^ We demonstrate 
the implementation of the method in realistic data sets in section]^ In section[^we present sensitivities to eccentric SMBBHs. 
In particular we compare two detection strategies towards the detection of eccentric binaries - a monochromatic search and 
a harmonic summing technique. Finally, we summarise and outline future work in section 


2 SINGLE-SOURCE GWS AND PULSAR TIMING 


In pulsar timing, the times of arrival (ToAs) of radio pulses from millisecond pulsars are measured and compared with 
predictions based on a timing model that describes the pulsar’s rotation (e.g., spin period and spin-down rate), the relative 
geometry between the pulsar and the observer, and other effects such as dispersion of radio waves due to interstellar medium 
( Lorimer fc Kramer]2005 Edwards et al.|2006 |. The differences between the measurements and predictions are called timing 
residuals. These residuals may be partly attributed to effects of GWs as they are not normally included in a timing model. 
Observations of a single pulsar cannot make a definite detection of a GW, but can lead to constraints on the strength of 


potential GWs (see, e.g., Yi et al. 2014 for a recent work). By timing an array of pulsars, GWs can be unambiguously 


detected by searching for correlated signatures among different pulsars. Typical PTA observations have a sampling interval 
of weeks and span over ~10 yr, implying a sensitive GW frequency range of ~1-100 nHz. 

In this work we are interested in single-source GWs, i.e., those coming from some particular directions in the sky. Timing 
residuals of such an origin can be generally written as: 

r{t, n) = F+{Cl)AA+(t) + Ex (n)AAx(t), (1) 

where U is a unit vector pointing from the GW source to the observer. The two functions F+{Q) and Fx (Cl) are the geometric 
factors as given by (e.g., |Lee et al. TMTt : 

1 


F+m = 


4(1 — cos 61) 


{(1-1- sin^ S) cos^ Sp cos[2(q: — Op)] — sin 25 sin25p cos(a — Op) -I- cos^ 5(2 — 3cos^ 5p)} 


( 2 ) 
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f’x(fi) = 


-{cos 5 sin 2Sp sin(a — Up) — sin 5 cos^ Sp sin[2(a — Op)]}, 


(3) 


2(1 — cosO) 

where cos 6 = cos5cos(5pCOs(Q — Qp) -I- sin <5 sin Jp with 9 being the opening angle between the GW source and pulsar with 
respect to the observer, a (ap) and 5 (5p) are the right ascension and declination of the GW source (pulsar) respectively. Note 
that we separate the polarization angle from F+ (fi) and Fx (Cl) and put it into AA+,x (f) for convenience of later analyses. The 
two functions F+{Cl) and Tx (12) are analagous to the antenna pattern functions used in the context of laser interferometric 
GW detection ( Thorne|1987| ). 

Because the wavelengths of GWs in the PTA band are much smaller than the pulsar-Earth distance, GW induced timing 
residuals are the combination of two terms - the Earth term A+,x(f) and the pulsar term A+,x(fp) (see, e.g., Jenet et al. 


20041: 


^^+,x (i) — ■A+,X (f) “ ^+,x (fp) 


(4) 


tp = t — dp{l — cos 9)/c, (5) 

where dp is the pulsar distance and we have adopted the plane wave approximation. A+{t) and Ax(t) are source-dependent 
functions. Throughout this paper, we only consider the correlated Earth-term signals. In fact for the case of continuous waves 
as expected from SMBBHs, pulsar terms will act as an extra source of uncorrelated noise for different pulsars. For GW bursts, 
whose duration is smaller than the data span, it is very unlikely that pulsar terms and Earth terms are simultaneously present 
in timing residuals for one particular pulsar unless the source sky direction is very close to that pulsar. Below we briefly 
discuss the signal properties for the three types of sources considered in this paper. 

At the leading Newtonian order SMBBHs in circular orbits are expected to generate pulsar timing signals with the 
following forms (see, e.g., Zhu et al. pOlit : 

ho 


A+{t) = [(1-1- cos^ l) cos 2‘ip sm(2n ft + (j>o) + 2 cos t sin 2ip cos{2n ft + ((iq)] 

210 f 

^x{t) = [(1 -b cos^ l) sin 2ip sm{2n ft + cj>o) — 2 cos l cos 2op cos(2io ft -b , 


( 6 ) 


(7) 


where the GW strain amplitude is ho = 2{GMcf'^^{iofY ^^with di, being the luminosity distance of the source and 
Me being the chirp mass defined as M®^^ = mim 2 (mi -b where mi and m 2 are the binary component masses, / is 

the GW frequency, i is the inclination angle of the binary orbit, tp is the GW polarization angle, (po is a phase constant. Note 
that we have neglected the frequency evolution over the observation span (typically ~ 10 yr). This represents the majority of 
circular binaries that are observable for PTAs as shown in Sesana & Vecchio ( |2010 1. 

Although it is well known that radiation of GWs circularizes the binary orbits |Peters fc Mathews 19631, the assumption of 
circular orbits is not always appropriate. Recent models for the SMBBH population including the effects of binary environments 
on orbital evolution suggest that eccentricity is important for GW frequencies < 10“® Hz ( Sesana|2013a R^^^etalj2014 1. In 
this work we use the expressions of GW-induced timing residuals for eccentric SMBBHs as given in Jenet et al. (20041. It is 


interesting to note that eccentric binaries emit GWs at multiple harmonics of the binary orbital frequency. At low eccentricities 
the emission is dominated by the second harmonic, while for high eccentricities the orbital frequency itself will dominate. 

Generally a GW burst is defined as a transient signal with a duration smaller than the observation span. This may be 
the only information we know about the source. For this reason we use a simple but general sine-Gaussian model for timing 
residuals induced by GW bursts: 

^2 ' 

( 8 ) 


A+{t) = A exp ^ {(1 -b cos^ t) cos2V’ cos[27r/o(t — to) + 4>o] — 2 cos t sin 2^/’ sin[27r/o(f — to) -b flio]} 

Ax{t) = A exp ^^ 2 °^ ^ {(1 -b cos^ t) sin2'!/)cos[27r/o(t — to) + epo] + ‘2> cos t cos 2^ sin[27r/o(t — to) -b ((> 0 ]} , 


(9) 


where A is the signal amplitude (in seconds), r is the Gaussian width, t is the source inclination angle, fo is the central 
frequency, to and tpo is the time and phase at the midpoint of the burst respectively. The sine-Gaussian model used here can 


represent qualitatively the signals for parabolic encounters of two massive black holes as studied in Finn & Lommen (20101 


and more recently in Deng (20141. The detection method that we will describe in the next section should apply equally well 


to other burst sources such as cosmic (super)string cusps and kinks (e.g.. 

Vilenkin||1981 Damour & Vilenkin||2000 Siemens 

et al.|2006l and triplets of supermassive black holes (|Amaro-Seoane et al. 

2010 b 


3 THE DATA ANALYSIS METHOD 

In this section, we describe how the singular value decomposition (SVD) can be used in a general method for the detection 
and localization of single-source GWs using PTAs. The method is adapted from the coherent network analysis method used 
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in the context of ground-based GW interferometers. There are some important features that are unique to PTA observations, 
e.g., (1) PTA data are irregularly sampled in contrast to nearly continuous sampling for ground-based experiments; and (2) a 
least-squares fitting process is performed to obtain estimates of timing parameters such as the pulsar’s spin period and its first 
time derivative, pulsar position and proper motion, etc. As we show later in this section, the SVD method proposed here relies 
on transforming the timing residuals of each pulsar to the frequency domain. For the idealized simulations (assuming even 
sampling and white Gaussian noise with equal error bars) used in section]^ a discrete Fourier Transform was used. In section 
for more realistic data sets we adopt a maximum-likelihood-based method to estimate Fourier components of the timing 
residuals making use of the noise covariance matrix (e.g., section 3 in Zhu et al.|2014 1. This is equivalent to the least-squares 
spectral analysis method (see, e.g., Goles et al.|2011 l. Our method works with post-fit timing residuals, i.e., after fitting ToAs 
for timing models of individual pulsars. The effects of the fitting process on our results will be discussed in section 

For a given source direction, the timing residuals from an array of Np pulsars can be generally written in the frequency 
domain as: 


dk = FkAk -f Hk, 


( 10 ) 


where the index k denotes the fe-th frequency bin, dk are timing residual data, Fk is the response matrix, Ak are GW 
waveforms and nk is the timing noise. The data are whitenecQso that 


( 11 ) 


where crfj, is the noise variance of the i-th pulsar at the fc-th frequency bin. Here dk, Ak and nk are all complex vectors, while 
the real whitened response matrix Fk is defined as 


dk = 

dlkjoik 

d2k/o2k 

> Ak = 

1 1 

+ X 

, nk = 

n-ife/cTifc 

n2fe/o'2fc 


b 

_1 




riNpk/o-Npk 


Fk = 


Fj^/nife 

/ark 

F2 

F}/a2k 

1- 

q 

F^J^Np 


( 12 ) 


For simplicity we will hereafter suppress the index k while keeping in mind that equations (|10p2 1 apply to each frequency 


bin in the analysis. Data for all frequencies can be stacked (e.g., in order of increasing frequency) to preserve the format 
of equation ( |10[ ), in which case F is a block-diagonal matrix (with a dimension of NpNf^ x 2Afc where Nk is the number of 
frequency bins) of Fk. 


The SVD of the response matrix F can be written as 


F = USV*,S = 


Si 0 

0 S2 

0 0 

0 0 


(13) 


where U and V are unitary matrices with dimensions of Np x Np and 2x2 respectively, and the symbol * denotes the 
conjugate transpose. Singular values in S are ranked such that Si ^ S 2 ^ 0. We can then construct new data streams as 
follows: 


d = U*d, A = V*A, n = U*n. 


(14) 


One can find that d = SA + n, and explicitly 


51 (V*A)^ -I- ni 

52 (V*A)2 + n2 

d= "3 


(15) 


nATp 


In the absence of GW signals, the real and imaginary parts for each element of d independently follow the standard Gaussian 


^ Here it is assumed that noise for different pulsars is uncorrelated. If not the full covariance matrix should be used in a whitening 
process. 
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distribution. Here di and d 2 are referred to the two signal stream^ since they contain all information about GW signals if 
present, while the remaining terms are ‘null streams’ (denoted by dnuii) since they have a null response to GWs. It has been 
shown that, in the context of GW detection with ground-based laser interferometers, null streams can be used as a consistency 
check on whether a candidate GW event is produced by detector noise or by a real GW l| Wen fc Schutz|2005 1 and ‘semi-null 
streams’ (i.e., d 2 if Si ^ S 2 ) can be included to improve the angular resolution ( Wen|2008 Wen et al.|2008 1. In this paper we 
only consider Earth terms in our signal model so the null streams (as constructed in the current way) are indeed ‘null’ but in 
reality they would have some response to the pulsar-term signals. In a future study pulsar terms will be incorporated in the 
detection framework with the addition of Np free parameters for pulsar distances in the response matrix F. 

It is straightforward to show that the maximum likelihood estimator for the GW waveform is 


A = Fd, where F = VSU*, S = 


1/si 

0 


0 

1/S2 


(16) 


The matrix F is the Moore-Penrose pseudoinverse of the response matrix F. The covariance matrix for the estimated waveform 
is 


r(A) = (^VS^SV*y 


(17) 


It is interesting to note that the statistical uncertainties of the estimated waveforms are a linear combination of and . 
Estimation of physical parameters can be obtained based on the estimated waveforms A and we leave this to a future study. 

Now we propose our detection statistics for the three types of signals discussed in the previous section. For monochromatic 
GWs, the detection statistic can be written as: 

Pmon = idl|"-f |d2|". (18) 

It is important to note that: (1) in the absence of GW signals, "Pmon follows a yf distribution with 4 degrees of freedom; and 
(2) the detection statistic applies to a given GW frequency and source sky location. In practice when such information is 
unknown, a search is usually performed to find the maximum statistic. 

For signals produced by eccentric binaries, we use a harmonic summing technique for which the detection statistic is 
given by: 

Pecc = ^ (ldl.,,J^ + id2.,fej^) , (19) 

where Nh corresponds to the highest harmonics considered in the search, fco is the bin number for the binary orbital frequency. 
In the absence of GW signals, Vecc follows a distribution with ANh degrees of freedom. In practice Nh should be determined 
as the one that gives the lowest false alarm probability (FAP). When the orbital frequency is unknown one should search over 
all possible frequencies to find the maximum statistic. 

For GW bursts with unknown waveforms, we adopt a time-frequency strategy in which the PTA data are first divided 
into small segments and then SVD is applied to each segment to output the two signal streams. Note that it is usually 
necessary to allow overlap between successive segments to catch signals occurring in the beginning or end of the segment. The 
detection of GW bursts of unknown waveforms generally involves searching for any ‘tracks’ or ‘clusters’ of excess power in 
the time-frequency space (e.g., Anderson fc Balasubramanian||1999 Wen &: Gair|20d5 |. So the detection statistic is designed 
as accumulating signal power in (a ‘box’ of) the time-frequency domain and can be written as: 

i/2 m/2 

■Pcwb)*, k) = E E (idb( i + a){k + b)\^ + |d2,(i + o)(fe + l,)|^) , (20) 

a= — l j1 b= — mf2 

for the i-th segment and fc-th frequency bin. Here I and m is the length of box in time and frequency respectively. If the data 
consist of only Gaussian noise, "Pcwb follows a yf distribution with ANh degrees of freedom (where Nh = I x m). 

The detection statistics proposed here all follow a noncentral distribution with their corresponding degrees of freedom 
when signals are present in the data. It is convenient to define the expected signal-to-noise ratio (p) as the noncentrality 
parameter 

(P) =iVdof+ p", (21) 

where the brackets (...) denote the ensemble average of the random noise process, Adof is the number of degrees of freedom for 
the ‘central’ distribution of V. We use p to quantify the signal strength in our simulations. Note that p^ equals the detection 
statistic calculated for noiseless data. In a frequentist detection framework, we are interested in the FAP of a measured 
detection statistic Vo, i.e., the probability that V exceeds the measured value for noise-only data. For the aforementioned 


^ Here the indices 1 and 2 refer to the first two terms in the new data streams d. More indices are used for time and frequency when 
necessary. 
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methods the single-trial FAP is given by 1 — CDF('Po;x^) where CDF(;x^) denotes the cumulative distribution function 
(CDF) for the central distribution in question. When a search is performed over unknown source parameters, the total 
FAP is 1 — [CDF('Pmax; for the maximum detection statistic Pmax found in the search with Wriai being the trials 

factor, which is defined as the number of independent cells in the searched parameter space for a grid-based search. 


3.1 Relation to the ‘A+Ax’ method 


In the pulsar timing software package TEMP02 ( [Hobbs et al.|[200^ , there exists a functionality with which two time series 
Al)?^ X (^) cs-n be estimated for a given sky direction. These two time series correspond to two polarizations of the coherent 
Earth-term timing residuals. Here we call it the ‘A+Ax’ method. In this method A+_x(f) are modelled as time-varying 
parameters and treated just like any other normal parameters in a timing model so that they can be estimated through a 
global least-squares htting routine. As observations for different pulsars are usually unevenly sampled and not at identical 
times, linear interpolation is used for such a global fit. The ‘A+Ax’ method enables one to simultaneously fit for single¬ 
source GWs and normal pulsar timing parameters. To avoid the covariance between the global fit for A+x (t) ^nd the ht for 
timing model parameters of individual pulsars, some constraints must be set on the two time series. Currently three kinds 
of constraints are implemented in TEMP02, namely, (1) quadratic constraints that correspond to pulsar spin parameters, 
(2) annual sinusoids for pulsar positions and proper motions, and (3) biannual sinusoids for pulsar parallax. These were first 
introduced and implemented in jKeith et al.| ( f2013| ) when correcting the dispersion measure variations for individual pulsars. 

The ‘A+Ax ’ method was hrst illustrated for arbitrary GW bursts in fig. 5 of |Hobb^ ( |2013| ). It was used in jZhu et al.| ( |2014| ) 
for an all-sky search for monochromatic GWs in the latest PPTA data set. A complete presentation on the ‘A+Ax’ method 
and its applications to single-source GW detection will be provided in a forthcoming paper (Madison et ah, in preparation). 
Here we briefly discuss the relation between the ‘A+Ax’ method and the SVD method proposed in this paper. Firstly the 


principle of both methods is identical since both use the response matrix F in a similar way and equation (161 essentially 


gives the least-squares solution to the two polarization waveforms, i.e., A in equation ( 161 is the frequency-domain equivalent 
of A+ X (^)' The critical difference is in the implementation - the ‘A+Ax ’ method works in the time domain whereas the SVD 
method works in the frequency domain. 

There are two advantages of the SVD method over the ‘A+Ax’ method: 

(1) The SVD method is much faster when searching over the whole sky is required. This is because it works with post-ht 
residuals after data for each pulsar have been fitted for a timing model and a search over unknown source sky location only 
involves doing the SVD for the response matrix F, while in the ‘A+Ax’ method a global ht is done for each searched sky 
location (as in the current implementation); 

(2) A truncated SVD may be used to improve the detection sensitivity, which applies to the case when si ^ S 2 , e.g., when the 
PTA has very low response to one of the two GW polarizations for some sky regions. This is possible especially for current 
PTAs whose sensitivities are dominated by a few best-timed pulsars. 

The ‘A+Ax’ method has been fully tested and implemented in real data for continuous GWs in Zhu et ah] ( [2014 1. We will 
demonstrate the effectiveness of the SVD method with some examples using idealized data sets in section [^ and then more 
realistic data sets in section 0 


4 EXAMPLES USING IDEALIZED DATA SETS 

For the purpose of illustration of our method, we consider an idealized PTA consisting of 20 PPTA pulsars. The simulated 
observations are evenly sampled once every two weeks with a time span of 10 yr. All observations contain stationary white 
Gaussian noise with a rms of 100 ns and equal error bars. The simulated data sets are produced as a combination of realizations 
of white Gaussian noise and Earth-term timing residuals. When we apply our detection statistics to noise-only data sets, it is 
confirmed that they follow a distribution with their respective numbers of degrees of freedom. 

In the following three subsections, we illustrate how our method works in terms of detection, localization and waveform 
estimation for the three types of signals considered in this work. The analysis is simplified as follows: (1) firstly the detection 
problem is demonstrated by evaluating the detection statistics at the injected source location; (2) then detection statistics 
are computed on a uniform grid of sky positions and for a range of frequencies to find the maximum; and (3) finally the 
frequency-domain waveforms A+(/) and Ax(/) are estimated at the actual source sky location. To show the correctness of 
the estimation process, a number of noise realizations were performed in order to obtain average estimates of A+,x (/) which 
were then compared with the true waveforms. For simplicity, in the relevant figures we only plot the absolute values of A+(/) 
and Ax (/), which are called spectral signatures. For the detection problem we calculate the single-trial FAP to quantify the 
detection signihcance, which is only appropriate when source parameters (such as sky location and frequency) are known as 
we assume for the following examples. The simulated signals are weak to moderately strong with signal-to-noise ratios ranging 
from 5 to 10. 
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Figure 1. Detection statistics (Pmon) as a function of frequency for a simulated data set that includes a strong monochromatic signal 
(solid black) and for the noise-only data set (dashed red). The vertical line marks the injected frequency (10 nHz), while the horizontal 
line corresponds to a single-trial FAP of 10“'*. 


Dec=90 



10 20 30 40 50 60 70 80 90 


Figure 2. Sky map of detection statistics calculated for a simulated data set that includes a strong monochromatic signal. The signal is 
simulated at the center of the map and the maximum detection statistic is found at “o” (which is the same as the actual source location). 
Sky locations of the 20 PPTA pulsars are marked with . 


4.1 Monochromatic waves 

For monochromatic waves, the simulated signal is characterized with the following parameters: ho = 1.16 x 10“^®, / = 10 nHz, 
cost = 1, Ip = 0, cpo = 0, (ct, S) = (0, 0). The expected signal-to-noise ratio is p = 10. Fig. shows the detection statistics 
(’Pmon) as a function of frequency evaluated at the injected sky location. The maximum statistic, which gives extremely strong 
evidence of a detection, is found at a frequency of 9.94 nHz. To localize this source, we calculate "Pmon for the same range of 
frequencies and for a grid of sky directions. At each sky direction, the maximum statistic over frequencies is recorded. Fig. 
shows that the source is successfully localized to where the signal is generated. 

Having already detected and localized the source, we use equations (16 ■ 171 to infer the GW waveforms. Fig. shows the 
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Figure 3. Estimated spectral signatures (thin solid black) for the same data set (white Gaussian noise + a monochromatic GW) used 
in Figs [T]and[2] along with true spectra (solid blue) and average estimates taken over 100 noise realizations (dashed red). 


estimated spectral signatures along with the true spectra, indicating a very good estimation as one would expect since the 
injected signal is very strong. To check the correctness of our method, we overplot in Fig. the average estimates taken over 
100 noise realizations - they are almost identical to the true spectra. 


4.2 Eccentric binaries 


The simulated source for this example is an eccentric SMBBH located in the Virgo cluster as described by the following 
parameters: mi = m 2 = 2 x 10 ®Mq, (11 = 16.5 Mpc, cost = 1, (q, S) = (3.2594, 0.2219), an orbital frequency of 5 nHz, an 
initial orbital phase of 1 and an eccentricity of 0.5. The signal is simulated using equations give in section 2 of |Jenet et al.| 
(20041. In order to ‘detect’ this simulated signal, we first apply our analysis at the injected sky location and experiment on 


the number of harmonics that should be considered. It turns out that a harmonic summation up to the second harmonic 
gives the lowest FAP. We will further discuss in section]^ the number of harmonics that should be included for binaries with 
different eccentricities. Fig. shows the detection statistics (Vecc) as a function of orbital frequency. The maximum statistic 
59.85, corresponding to a FAP of 5 x 10”'^°, is found at 4.87 nHz. The expected signal-to-noise ratio for this injection is p = 7 
when we perform a harmonic summing up to the second harmonic. Then we show in Fig. the detection statistics calculated 
for the whole sky (after maximizing over orbital frequencies for each sky direction). The maximum statistic is found at a grid 
point close to the injected sky location. 

Fig. [6] shows the estimated spectral signatures assuming that we know the actual source location. Since the detection 
statistic varies only slightly within an area of tens of square degrees (as shown in Fig.[^, similar results should be obtained if 
we apply the spectral estimation analysis in the sky location where the maximum statistic is found. It is shown that reasonably 
good estimates are obtained only for the first two harmonics. 


4.3 GW bursts 

In this example the simulated GW burst signal is described by the following parameters: A — 100 ns, r = 60 days, cos b = 0.5, 
/o=50 nHz, Ip = 1.57, cpo = 0, (a, 5) = (0.9512, —0.6187), i.e., originating from the Fornax cluster. The signal occurs in the 
middle of our observations (MJD 55250). Waveforms of A+,x(t) are shown in the upper panel of Fig.[^ For this burst source 
the amplitudes of timing residuals are < 100 ns, which means the signal should not be visibly apparent in individual pulsar 
data set. 
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Figure 4. Detection statistics {'Pecc) eis a function of orbital frequency for a simulated data set that includes a moderately strong signal 
(solid black) and for the noise-only data set (dashed red). The signal is produced by an eccentric SMBBH with an eccentricity of 0.5. 
Detection statistics are calculated as a harmonic summation up to the second harmonic. The vertical line marks the injected orbital 
frequency (5 nHz), while the horizontal line corresponds to a single-trial FAP of 10“^. 
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Figure 5. Sky map of detection statistics calculated for a simulated data set that includes a simulated signal produced by an eccentric 
SMBBH located in the Virgo cluster (marked by a white diamond). The maximum statistic is found at “o”. Sky locations of the 20 
PPTA pulsars are marked with . 


To dig out this burst signal from noisy data, we divide the 10-yr data set into segments of length of 300 days and 
calculate the detection statistic given by eqnation 20 in the time-frequency domain. The segment length is (roughly) chosen 
based on the knowledge of /o, i.e., having > one cycle per segment. In practice since the time-frequency analysis for PTA 
data shonld not be limited by computational power, many trials of segment length can be performed to search for signals of 
different durations. The lower panel of Fig.j^shows results of the time-frequency analysis assuming known source sky location. 
The most signihcant statistic 38.36, corresponding to a FAP of 9 x 10“®, is found at MJD 55257 and a frequency of 51.67 
nHz. Therefore our analysis has clearly detected the injected signal and correctly identified its occurrence time and central 
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Figure 6. Estimated spectral signatures (thin solid black) for the same data set (white Gaussian noise + a signal expected from an 
eccentric SMBBH) used in Figs and along with true spectra (solid blue) and average estimates taken over 100 noise realizations 
(dashed red). 


frequency. For this example, the signal is very ‘isolated’ in the time-frequency space, so it is sufficient to only look at the 
maximum statistic in the time-frequency map (in this case, the signal-to-noise ratio is p = 5). 

Fig .[^ shows the detection statistics evaluated at the whole sky (after maximizing over frequencies for each sky direction). 
The maximum statistic is found at a grid point close to the injected source location. Fig. shows the inferred spectral 
signatures compared against the true spectra. As the injected signal is relatively weak for this example, the spectrum is not 
recovered as well as the previous two examples. For both plots only the central segment that contains the majority of signal 
power (as illustrated in Fig. was used. 


5 MORE REALISTIC DATA SETS 


Examples given in the previous section have assumed idealized PTA observations that are evenly samplecBand contain only 
stationary white Gaussian noise with equal error bars. Realistic features typical to real PTA data include: (1) observations 
are irregularly sampled and have varying ToA error bars. It is also common that the data span varies significantly among 
different pulsars; and (2) low-frequency (‘red’) timing noise may be present for some pulsars. To simulate realistic data sets, 
we make use of the actual data spans, sampling and error bars of the PPTA 6-yr Data Release 1 (DRl) data set that was 


published in Manchester et al. (20131 and used in Zhu et al. (20141 for an all-sky search for continuous GWs. The PPTA DRl 
data set is publicly available online at a permanent linl^ As a check, in some simulations an uncorrelated red-noise process 
with a power-law spectrum is also included. This has no effect on the results since we assume the noise spectrum is known 
and thus include it explicitly in the noise covariance matrix. In actual analyses of real data, the noise estimation is a very 
important step and we do not attempt to address it in this work. 

Fig. shows the distribution of the detection statistics "Pmon i n th e presence of signal along with the distribution of 
corresponding squared null streams (|dnuiip) as defined in equation (15l. The injected signal is characteristic of a circular 
SMBBH described by the following parameters: ho = 1 x 10“^^, f = 20 nHz, cost = 0.6172, ip = 4.1991, (pa = 1.9756, 


® Simultaneous to the estimation and correction of the dispersion measure variations for multiple-band pulsar timing data ||Keith et al.| 
one obtains estimation of the common-mode signal that is independent of radio wavelengths. These common-mode data sets are 
evenly sampled and can be searched for GWs. 
http: //dx. dol. org/10.4225/08/534CC21379C12 
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Figure 7. Upper panel; simulated timing signals from a GW burst. Lower panel: detection statistics Powb (defined in equation |20[ | 
calculated in the time-frequency domain. 
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Figure 8. Sky map of detection statistics calculated for a simulated data set that includes a simulated GW burst originating from the 
Fornax cluster (marked by a white diamond). The maximum statistic is found at “o”. Sky locations of the 20 PPTA pulsars are marked 
with . 


4.1 


To obtain 


(a, 5) = (0, 0). The expected signal-to-noise ratio is p = 10, which is identical to the example shown in section 
the distribution of the detection statistic and null streams, we perform 10"^ realizations of white Gaussian noise and for each 
noise realization keep the search parameters (i.e., source frequency and sky direction) fixed at their injected values. We can 
see that both match the expected distributions perfectly well. As mentioned in the previous section, null streams can be used 
as a consistency check in the case of a detected GW signal - if the detected signal is due to a GW, null streams should follow 
the expected noise distribution while otherwise false alarms caused by any noise processes are very unlikely to exhibit such a 
property. 

Fig.f^shows the detection statistics as a function of frequency calculated for realistic simulated data sets in the absence 
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Frequency (FIz) 


Figure 9. Estimated spectral signatures (thin solid black) for the same data set (white Gaussian noise + a GW burst) used in Figs[^ 
and|^ along with true spectra (solid blue) and average estimates taken over 100 noise realizations (dashed red). 


or presence of the same monochromatic signal as in Fig. |10[ We have considered two statistics - "Pmon proposed in this paper 
and the J^e-statistic that was hrst proposed in Babak & Sesana (20121 and later strengthened in Ellis et al. (20121, both 
giving nearly identical results. In the signal-present case, the maximum value (129.4) of "Pmon is found at 20.5 nHz, while 
the largest J^e-statistic (125.4) is measured at 20.6 nHz. It is also obvious that at some high frequencies there is a spectral 
leakage problem that applies to both methods. This is not surprising given the similar principle of the two statistics: (1) the 
J^e-statistic works fully in the time domain, but it also involves a process of maximum-likelihood estimation of the fourier 
components in individual pulsar data set; (2) in the derivation of the J^e-statistic, the four (time-dependent) basis functions 
of the monochromatic timing residuals are essentially equivalent to the two orthogonal sine-cosine pairs of j4-i-,x (!)• 

Fig. shows the results of an all-sky search using the SVD method for the same signal-present data set as used in Fig. 
[m The maximum detection statistic 130.1 is found at (a, 3) = (6.185, —0.133) with a frequency of 20.5 nHz. The source is 
localized to a direction close to the injected location. However, compared with Fig. we can see that the angular resolution 
for realistic PTAs is much worse because the sensitivity is dominated by a few good ‘timers’ in the array. It is worth pointing 
out that the SVD method is also advantageous in terms of computational efficiency over that of the 7^e-statistic, since for the 
latter the amount of computation is proportional to the total number of sky locations being searched. 

Similar results were obtained if the simulated data sets used in Figs 1 1 0f[T^ have gone through the TEMP02 fitting process 
for a full timing model for each pulsar. In fact, it has been shown that the fitting process can be approximated by multiplying 
the timing residuals by a data-independent and non-invertible linear operator matrix (e.g., Demorest|2007 1. This matrix can 
be explicitly included in the noise covariance matrix and does not pose a problem for our analysis method. The timing model 
ht is known to absorb some GW power at low frequencies (close to l/Tspan with Tspan being the data span, because of the 
ht for pulsar spin period and its first time derivative), and in some narrow bands centred around 1 yr“^ (the fit for pulsar 
positions and proper motions), 2 yr”'^ (the fit for pulsar parallax) and high frequencies that correspond to binary periods for 
pulsars in binaries (e.g., hg. 1 in Cutler et al.|2014). These effects are similar to those caused by applying constraints on the 


X (f) time series as discussed in section 


3.1 


6 SENSITIVITIES TO ECCENTRIC SMBBHS 

In this section we are interested in current PTA sensitivities to eccentric binaries. Given that previous published PTA searches 
for GWs from SMBBHs have assumed zero eccentricity, we attempt to answer the following question - when is a monochromatic 
search sufficient to detect eccentric binaries? For this purpose, we use simulated data sets that have the same sampling and 
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(a) 



Figure 10. (a) Probability distribution of the detection statistic T^mon (solid blue) in the presence of a monochromatic signal compared 
to the expected noncentral distribution (dashed red) with 4 degrees of freedom and an noncentrality parameter = 100. The 
simulation was performed with search parameters fixed at the injected values and 10^ realizations of white Gaussian noise, (b) For the 
same simulation as panel (a) but instead showing the probability distribution of the squared null streams |dnuiiP (solid blue) compared 
to the expected distribution with 2 degrees of freedom (dashed red). 
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Figure 11. Detection statistics as a function of frequency for realistic simulated data sets in the absence (thin curves) and presence 
(thick curves) of a strong monochromatic signal. Here we compare two statistics - "Pmon (solid red) proposed in this paper and the 
J^e-statistic (black dash). All statistics are evaluated at the injected source sky direction. The vertical line marks the injected frequency 
(20 nHz), while the horizontal line corresponds to a single-trial FAP of 10“"^. 


Dec=90 



20 40 60 80 100 120 


Figure 12. Sky map of detection statistics calculated for a simulated realistic data set that includes a strong monochromatic signal. 
The source is simulated at the center of the map (indicated by a white diamond) and the maximum detection statistic is found at “o”. 
Sky locations of the 20 PPTA pulsars are marked with . 


error bars as the PPTA DRl data set, and consider two detection statistics "Pmon and Vecc- The sensitivities are parameterized 
by the luminosity distance within which 95% of sources are detectable at a fixed FAP of 10“'*. 

To facilitate our calculations, we simulate signals due to eccentric binaries drawn from uniform distribution in cos 5 and 
a and evaluate both "Pmon and Vecc for noiseless data at the injected sky location. For the demonstration here, we consider 
sources with a chirp mass of 10 ®Mq and an orbital frequency of 10 nHz with other parameters such as cos t, initial orbital 
phase and polarization angle randomized. Since here the data consist of signals only, the detection statistic equals the squared 
signal-to-noise ratio p® and scales inversely with d® . We use this scaling relation to find the value of di that corresponds 
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Figure 13. Luminosity distance (di), within which 95% of eccentric SMBBHs with a chirp mass of 10 ®Mq and an orbital frequency 
of 10 nHz could be detectable with the current PTA sensitivity, as a function of orbital eccentricity (eg). Two detection strategies 
are considered: a monochromatic search (black dots) and a harmonic summation technique (blue open circles). Note that the sensitive 
distance scales as 


to the given detection threshold. We perform 10^ Monte Carlo simulations and choose the 95% quantile as a point in the 
sensitivity plot. 

Fig. shows the sensitivity to GWs from eccentric SMBBHs using two detection strategies - a monochromatic search 
and a harmonic summing technique. It is clearly demonstrated that for high and low eccentricities a monochromatic search 
is better than the harmonic-summing search, while the latter is more sensitive in the moderate regime (0.3 < eg < 0.55). 
This is in good agreement with the fact that for high and low eccentricities the GW emission is dominated by the orbital 
frequency and its second harmonic respectively. In the high- and low-eccentricity regimes, adding incoherently power from 
secondary harmonics is not beneficial because when a harmonic summing technique is adopted the detection threshold should 
be increased for a fixed FAP. We also find that the inclusion up to the third harmonic is sufficient for all possible eccentricities 
in this example. 


7 SUMMARY AND FUTURE WORK 


PTA experiments have steadily improved their sensitivities over the past few years and produced very stringent constraints on 
the evolution of ensembles of SMBBHs. While the traditional focus was a stochastic background from SMBBHs throughout 
the Universe, possibilities of detecting and studying single-source GWs using PTAs have been explored in depth in recent 
years. Although it may be more challenging, detections of individual sources will provide rich information on the properties 
of the GW emitters such as the orbital eccentricities, masses and even spins of SMBBHs. This prospect becomes increasingly 
important as the Chinese 500-meter aperture spherical telescope (FAST) is expected to come online in 2016 (Nan et al.|2011 


[Hobbs et al.||20l4| | and the planned Square Kilometer Array (SKA) will be operational within a decade ( [Lazio [[2013[ ). Both 
FAST and SKA will provide a major step forward in terms of single-source GW detection with PTAs. 

Based on the coherent network analysis method used in ground-based interferometers, we proposed a method that is 
capable of doing detection, localization and waveform estimation for various types of single-source GWs. We demonstrated 
the effectiveness of this method with proof-of-principle examples for GWs from SMBBHs in circular or eccentric orbits and 
GW bursts. We also demonstrated the implementation of this technique using realistic data sets that include effects such 
as uneven sampling and varying data spans and error bars. The new method is found to have the following features: (1) 
it is fast to run especially for all-sky blind searches; (2) it performs as well as published time-domain methods for realistic 
data sets; and (3) null streams can be constructed as a consistency check in the case of detected GW signals. Finally, we 
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presented sensitivities to eccentric SMBBHs and found that (1) a monochromatic search that is designed for circular binaries 
can efficiently detect SMBBHs with both high and low eccentricities; and (2) a harmonic summing technique provides better 
detection sensitivities for moderate eccentricities. 

Our future work will be: 

(1) to develop a fully functional data analysis pipeline that can be tested in future IPTA mock data challenge^ and applied 
to real data; 

(2) to make comprehensive comparisons (in terms of detection sensitivity, parameter estimation and speed) between the 
method proposed in this paper and published methods such as the ‘A+Ax’ method and the J^e-statistic. Again this could be 
done along with a future IPTA data challenge; 

(3) to include pulsar terms in the detection framework for continuous GWs and investigate how the information of pulsar 
terms can be exploited to improve the PTA angular resolution. 
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